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Abstract 

Performance of Wang-Landau (W-L) algorithm in two continuous spin models is 
tested by determining the fluctuations in energy histogram. Finite size scaling is 
performed on a modified XY-model using different W-L sampling schemes. Difficul- 
ties faced in simulating relatively large continuous systems using W-L algorithm are 
discussed. 
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The Wang-Landau (W-L) algorithm for Monte Carlo (MC) simulation, introduced in 
2001 [Ij has since been applied to a wide range of problems in statistical physics [21 [3l HJ O [6] . 
In most of these investigations the authors have applied the W-L algorithm to systems 
with discrete energy levels which include the Ising and the Potts model. Relatively fewer 
papers have so far appeared on continuous models [H [HI [9] and in such systems one uses 
a discretization scheme to divide the energy range of interest into a number of bins which 
label the macrostates of the system. 
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In the W-L algorithm, one directly determines the density of states fi(-Ej) (in practice, 
its logarithm g{Ei)) of the system with i=l,2,..n being the bin index. These macrostates 
are sampled with a probability which is inversely proportional to the current value of the 
density of states. If a trial move passes the probability test, g{Ej) for the new energy 
Ej is modified as g{Ej) — > g{Ej) + In/, where the modification factor / is > 1 in the 
beginning of the simulation. In case a trial move fails, the density of states corresponding 
to the old value Ei of the energy is modified. A histogram record H{Ei) of all states 
visited is maintained throughout the simulation. When the g{Ei) corresponding to a certain 
macrostate is modified by adding In / to it, the corresponding H{Ei) is modified as H{Ei) 
H{Ei) + 1. In the beginning of the simulation the g{Ei)'s for all macrostates are initialized 
to zero. In the original version of W-L algorithm, an iteration is considered to be complete 
when the histogram satisfies a "flatness" criterion. This means that H{Ei) for all values of i 
, has attained 90% (or some other preset value) of the average of all H{Ei). In the following 
iteration / is reduced, the H{Ei)^s are reset to zero, and the process is continued till In/ 
is as small as 10^^ or 10^^. Since the history of the entire sampling process determines 
the density of states, the W-L algorithm is non-Markovian besides being multicanonical in 
nature. 

In course of the random walk in a W-L simulation, the fluctuation of the energy his- 
togram, for a given modification factor /, initially grows with the number of Monte Carlo 
sweeps and then saturates to a certain value. Because of the nature of the W-L algorithm, 
which has been described above, the value of the histogram fluctuation determines the error 
which is generated in the resulting density of states. Zhou and Bhatt (TU] carried out a 
mathematical analysis of the W-L algorithm. They proved the convergence of the iterative 
procedure and have shown that the error in the density of states, for a given /, is of the 
order of -y/(ln/). This finding has been tested by Lee et.al [H] who performed extensive 
numerical test in two discrete models. In Monte Carlo simulation on two two-dimensional 
Ising models, namely, the ferromagnetic Ising model (FMIM) and the fully frustrated Ising 
model (FFIM), Lee et.al have shown that the fluctuation in the histogram increases during 
an initial accumulation stage and then saturates to a value which is inversely proportional 
to y^ln/) and they were of the view that this feature is generic to the W-L algorithm. As 
is shown in references [10] and [11] and also in the later part of this Letter, the resulting 



error in the density of states is then of the order of y {In f), which is in agreement with the 
prediction of Zhou and Bhatt. 

Methods which are alternatives to the requirement of the histogram flatness, for decid- 
ing where to stop an iteration for a given modification factor, have been proposed in the 
work described in references [10] and [H]. Zhou and Bhatt [lOj were of the opinion that 
an iteration may be stopped when the minimum number of visits to each macrostate is 
l/^\nf). On the other hand, Lee et.al [H] proposed that an iteration can be stopped 
when the number of Monte Carlo sweeps, for a given value of /, is such that the satura- 
tion of the histogram fluctuation has been reached, since continuing the simulation for this 
particular modification factor is unlikely to reduce the error in the density of states any 
further. Subsequently, there have been a number of proposed improvements and studies of 
the efficiency and convergence of this algorithm [121 [13], [M] . 

In this Letter we first present the results of our investigation on the growth of histogram 
fiuctuations in two continuous models. One of our aims is to check if the conjecture of 
Lee et.al., that the nature of the dependence of the maximum of the histogram fiuctuation 
on the modification factor / is model independent, can be extended to lattice spin models 
with continuous energy spectrum. For this purpose, we have chosen a two dimensional spin 
system, where the spins, confined to a square lattice and free to rotate in a plane, say the 
xy- plane, (having no z-component) interact with the nearest neighbours via a potential, 



where 6 is the angle between the interacting spins and is a parameter to be chosen. This 
model, now known as the modified XY- model, was first introduced by Domany et.al in 
1984 [15j. For p"^ = 1, the model is simply the conventional two-dimensional XY- model, 
which is known to exhibit a quasi-long-range-order- disorder phase transition mediated by 
the unbinding of topological defects |T6]. For large values of p^, the potential of eqn.(l) 
has a sharp well structure for small values of 6 and the model exhibits a rather strong first 
order phase transition. Thus, the model with interaction determined by eqn.(l), for values 
of = 1 and 50, effectively behaves like two different models, having characteristically 
different phase transitions, although the symmetry of the Hamiltonian and the lattice and 
the spin dimensionalities remain the same. A similar change in the nature of the phase 
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transition has been observed in the two dimensional Lebwohl-Lasher (LL) model and a 
modified version of it [T7j. In this model, with spin dimensionality n = 3 and lattice 
dimensionality d = 2, the nearest neighbour spins interact via a potential —P2{cos6ij) and 
in the modified (LL) model this is —P4{cos6ij) where, P2 and P4 are the second and the 
fourth Legendre polynomials respectively. Both models have the 0(3) as well as the local 
Z2 symmetry. But as the nature of the interacting potential is modified, the transition 
changes from continuous to a sharply first order one. It has also been observed that for the 
n = 3 and d = 3 LL model, if one adds a P4 interaction to the usual P2 term, the nature of 
the phase transition changes from a weakly first order one to one where the discontinuities, 
which characterize the first order transition, grow sharply [HI [191 1201 [211 [22l [23] . 

The p"^ = 50 model has very large fiuctuations in energy, which is manifested in a huge 
peak in the specific heat. This is accompanied by very long relaxation times that make it 
difficult to obtain good statistics in the simulation using the conventional canonical sampling 
as is done in the Metropolis method [24] • In their original work, Domany et.al [15] had to 
carry out the simulation of this model in the roughening representation [25] and employ 
long runs. 

Using different approaches mentioned below, we have also tested the performance of the 
W-L algorithm in simulating the = 50 model. In the first method, which we call 'a', 
a pre-assigned number of MC Sweeps, determined from the saturation of the histogram 
fiuctuation, is used to stop an iteration with a given modification factor. In the other 
two methods, 'b' and 'c', a minimum number of visits to each macrostate (bins) namely 
l/y^ln/) and l/(ln/) respectively have been used. So, in our work three different criteria, 
instead of testing the flatness of the histogram, have been used for stopping a particular 
iteration. Besides the methods 'a' and 'b' suggested in references [TT] and [10] respectively, 
we have chosen method 'c' because the proposition of l/-y/(ln/) visits to each macrostate 
does not directly follow from the fact that the error in the density of states is of the order of 



Y ( In /). We have exercised this option, although it leads to a huge increase in the amount of 
computation, to see in particular if it leads to any improvement in the results of simulation. 
This is necessary particularly in view of the fact that, we have worked in a continuous model 
and the number of microstates which correspond to each bin is enormously large. 

Before presenting our results we explain in detail the notations and symbols relevant to 




the problem. We represent by Pk, the saturation value of the energy histogram fluctuation in 
the k^^ iteration. Let fk be the modiflcation factor for the k^^ iteration. One usually starts 
with a modiflcation factor / = /i > 1 and uses a sequence of decreasing /^'s (k=l,2,3,....) 
deflned in some manner. One MC sweep is taken to be completed when the number of 
attempted single particle moves equals the number of particles in the system. The error in 
the density of states after the n*^ iteration has been performed, is directly related to /3j for 
i > n, the saturation values of the fluctuations. 

In the W-L algorithm the logarithm of the density of states after n iterations is given by 

n 

g^{E,) = J2H,{E,) ln(/fc) (2) 

k=l 

where, H^^Ei) is the accumulated histogram count for the i^^ energy bin during the k*^ 
iteration. In order to get an idea of the fluctuations in the histogram and its growth with 
the number of MC sweeps we subtract the minimum of the histogram count Wf, which occurs 
in the histogram after the j*^ MC sweeps has been completed during the k^^ iteration, i.e 
we consider the quantity 

Hm) = Hm)-hi (3) 

(It may be noted that hi. does not refer to any particular bin and may occur randomly in 
any of the bins). 

The quantity Hl{Ei) is now summed over all bins to give ^H^. 

AHi=j:Hm) (4) 

i 

AHl is thus a measure of the fluctuations which occurs in the j^^ MC sweeps during k^'^ 
iteration and is a sort of average over all macrostates or bins. AiJ^ fluctuates with j 
because of statistical errors and its mean value taken over j is nothing but (3k- The error 
in the logarithm of the density of states, summed over all energy levels or bins, after the 
completion of n iterations is therefore given by 

oo 

Vn= J2 ln(/fc) (5) 

k=n+l 



Since the predicted value of the error rin is of the order of y ( In /„) [TU], one expects that the 
histogram saturation value for the n*'' iteration, should be proportional to l/i/(ln/„). 



In the present work we have determined for the two lattice models we have defined, the 
dependence of the quantity Ai7^, given by eqn(4), on j, the number of MC sweeps for a 
given iteration denoted by k. Besides this, by using lattices of linear dimension L = 8 to 80 
for the — 50 model we have performed finite size scahng (FSS) using methods 'b' and 'c' 
so that the relative performance of the two methods can be judged. 

For both the models and for the purpose of testing the fiuctuation in histograms we 
have worked on two system sizes, namely 8x8 and 16 x 16, and nearest neighbour inter- 
actions along with periodic boundary condition were always used. The starting value of 
the modification factor In/i was taken to be 0.1 and the sequence ln/„_|_i = ln/„/(10)^/^ 
was chosen and for the purpose of determination of fluctuations the minimum In / used was 
10~^. Clearly, the chosen sequence of / is to ensure that it gets reduced by a factor of 10 
after four iterations. We have determined the quantity Ai7^ defined by eqn(4) at intervals 
of 10^ MC sweeps and the maximum number of sweeps chosen for a given value of / is such 
that the saturation of the histogram is clearly evident. The energy range accessible to the 
systems we have investigated is to 4 per spin. To reduce the computer time we have taken 
the system energy range to be 3 to 253 for the 8x8 lattice and 10 to 1004 for the 16 x 16 
one. The width of the energy bins was taken to be 0.2 in all cases. 

In fig. 1, we have plotted the fiuctuations in the histogram Aif^ against the number of 
MC sweeps, j for four values of the refining factor /. The plots shown are for the 16 x 16 
lattice for = 50 and for In / equal to 10^^, 10^'^, 10^"^, 10~^. We did not go to values 
of In/ less than 10~^ as it takes a very large CPU time. Averages were taken over fifty 
independent simulations to improve the statistics. Similar plots were also taken for the 
p'^ — 1 model. It is evident from figure 1 that Aif^ increases initially and then saturates 
and as / gets smaller, the saturation value as well as the number of MC sweeps necessary 
to reach the saturation increases. The standard errors calculated from the fifty independent 
simulations are also shown. In fig. 2 we have plotted the logarithm of the saturation value, 
Pk ^-e ln( Pk) vs la{lnf) for the p'^ — 1 model for 8 x 8 and 16 x 16 system sizes and those 
for the = 50 model are depicted in fig. 3. The error bars are not shown since they are 
of size smaller than the dimension of the symbols used for plotting. From these figures it is 
clear that 

/9feOc(ln/)" (6) 



where the index a = -0.508 ± 0.002 and -0.507 ± 0.002 for 8 x 8 and 16 x 16 respectively 
for ^2 = 1 and a = -0.500 ±0.003 and -0.504 ±0.003 for 8 x 8 and 16 x 16 respectively for 
= 50. This is in agreement with the prediction of Zhou and Bhatt [10] and the findings 
of Lee et.al [TT] for the two discrete Ising models. 

We denote by A the number of MC sweeps necessary to reach the saturation value of 
the histogram fluctuation. Our results show that In A is a linear function of ln(ln/), like 
InA = a±61n(ln/) 

i.e., X = exp{a){lnff (7) 

where a, b are parameters. We have plotted In A against ln(ln /) in figure 4 for the = 50 
for two lattice sizes. The fits yield a = 8.78 ± 0.18 and b = -0.526 ± 0.02 for the 8 x 8 
lattice. For the 16 x 16 system these parameters are a = 8.03 ± 0.33 and b = —0.580 ± 0.04 
respectively. It therefore seems that A depends on the system size rather weakly. 

During the W-L simulation as the modification factor becomes smaller, the changes in 
the density of the state profile also becomes finer. We denote by Agk{Ei) the change in the 
logarithm of the density of states resulting from the k^^ iteration and Ei is the bin energy. 
The plots of Agk{Ei) against Ei for four iterations in the 8 x 8, = 50 model is shown in 
fig. 5. As expected, the fluctuation in Ag decreases with the decrease in /. These diagrams 
reveal a little more information than that is evident from the AHl curves in fig. 1, as the 
fluctuation spectrum is here available as a function of the energy bin number, rather than 
the fluctuation summed over the bins. This is an important factor in every simulation since 
AHl '^^^ saturate at a given value in spite of not being distributed uniformly over all bins. 
That this does not happen for this system, particularly as / gets reduced, is apparent from 
these diagrams. However we have noticed that for bigger systems this behaviour of Agk{Ei) 
is violated. 

Before presenting results of FSS in the p"^ = 50 model, we outline the approach of 
Lee and Kosterlitz EZl who had worked out a method for doing FSS in systems with 
a first order phase transition. When simulating an unknown system, two issues need to 
be resolved. One is the nature of phase transition and the other is the determination of 
various thermodynamic quantities. Lee and Kosterlitz proposed a convenient method for the 
determination of the order of the phase transition and this can be applied to MC simulations 
in systems having linear dimension less than the correlation length. For a temperature driven 



first order transition in a finite system of volume L*^ with periodic boundary condition 
one needs to compute a quantity N{E; (3, L) wliicli is tlie liistogram count of tlie energy 
distribution. Tlie = 50 model has a characteristic double peak structure for N{E; f3, L) 
in the neighbourhood of the transition temperature. The two peaks of at Ei{L) and 
E2{L) corresponding respectively to the ordered and disordered phases are separated by a 
minimum at Em{L). A free energy like quantity is defined as 

A{E- (3, L,Af) = - In N{E; (3, L) (8) 

where M is the number of configurations generated. The quantity A{E\ (3, L, A/") differs from 
the free energy F{E] j3, L) by a temperature and A/" dependent additive quantity. A bulk 
free energy barrier can therefore be defined as 

AF{L) = A{E^^;[3,L,M)-A{Ei;[3,L,Af) (9) 

It may be noted that at the transition temperature, A{Ei; P, L,Af) = A{E2; P, L,J\f) and 
AF is independent of A/". For a continuous transition AF{L) is independent of L and for 
a temperature driven first order transition it is an increasing function of L, even when L is 
smaller than the correlation length, ^, prevailing at the system at the transition temperature. 
If one is in a region where L is much greater than ^, then AF obeys the scaling relation [28] 

AF ~ L"^-^ (10) 

For FSS analysis in the = 50 model we have used systems of size Lx L with the linear 
dimension L = 16, 32, 48, 64 and 80. We observe that while using the method 'a' for the W-L 
algorithm, the termination of the simulation for a given / at the pre-determined value of 
MC sweeps, as is necessary to reach the saturation value of the histogram fluctuation, does 
not work for L > 32. This is because, a large number of bins, in the low energy range are not 
visited at all. The problem becomes more severe as the system size increases. The error in 
the density of states resulting from the energy range not being sampled at all leads to wrong 
values of thermodynamic variables. This means that, mere specification of a pre- assigned 
value of the necessary number of MC sweeps, without monitoring if each bin is adequately 
sampled, does not provide us with a method which is an alternative to the flatness test of 
the histograms. For clarity we present in table 1, the energy ranges for different system 



sizes, which are not sampled at all for MC sweep ranging upto the estimated pre-assigned 
number. 

Next we turn to results of FSS using the methods 'b' and 'c'. Figure 6 shows the specific 
heat per particle plotted against the dimensionless temperature for L ranging from 8 to 
80 and the sharpness of the C^-peaks grow rapidly with increasing system size. The free 
energy A is plotted against energy/particle for all lattice sizes in figure 7. The growth of 
a free energy barrier between two wells can be seen. The data for both fig. 6 and 7 have 
been obtained by using method method 'b'. In fig. 8 we show the finite size scaling of the 
peak height of obtained by the two methods where has been plotted against L^. In 
accordance with first order scaling laws [28], the behaviour of Cy vs is linear and the 
difference between the results obtained from the two methods of sampling is insignificant. 
In figure 9, we have plotted the free energy barrier AF at the transition temperatures vs 
the lattice size L. Using method 'b' we get a good linear fit (denoted by filled circles in the 
figure) which is expected from the analysis of Lee and Kosterlitz. However, for the method 
'c' the linear behaviour of AF vs L is not obtained (denoted by filled triangles in the figure) 
although AF remains a monotonically increasing function of L, which also characterizes a 
first order phase transition [26l [27] . 

In figure 10 we present the finite size scaling behaviour of the transition temperatures 
for different values of L obtained from i) peak position of and ii) fine tuning of free 
energy with temperature till two equally deep minima are obtained. We find that the fit 
of the transition temperature vs is linear, thus obeying the finite size scaling rules for 
a first order transition. When we fit the data, we get two sets of straight lines for the two 
methods we used. Table 2 shows the thermodynamic limits of the transition temperatures 
thus obtained and the difference between the two methods appears to be insignificant. 

We end this letter by noting down our conclusion. In the first part of the work we have 
found that, in the two planar spin models with continuous energy spectrum, the fluctua- 
tion in the energy histogram, after an initial increase, saturates to a value proportional to 
l/-y/(ln/) where / is the modification factor in the W-L sampling. Since the error in the 
resulting density of states depends on this saturation value, it may seem that W-L sampling 
should be carried out only till the saturation value of histogram fluctuation is attained. We 
have shown that this does not work even for systems of moderate size, as an energy range 



near the ground state is not sampled at all. In the second part of our work, we have carried 
out the W-L simulation of the — 50 model for lattices of size upto 80 x 80 using two 
schemes, where the minimum number of visits to each bin are 1/^lnf) and 1/ln/ respec- 
tively. The second scheme, although requiring a vast amount of computation, particularly 
as / gets smaller, docs not give results which arc significantly different from the method 
employing l/^\nf) visits. All the expected FSS behaviour for first order phase transition 
has been seen to be obeyed in the simulations we have done and the l/y^ln/) method 
seems to be adequate. It may be noted that even the l/y^ln/) visit scheme in W-L sim- 
ulation in lattices of bigger size seems to be impractical in view of the huge computer time 
necessary. This is a serious problem with the W-L algorithm when applied to continuous 
lattice spin models. 
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Table 1: The energy ranges used in the simulation for the — 50 model for system sizes 
of linear dimension L using the method of pre-assigned MC sweeps (method 'a'). The third 
column shows the energy range which is not visited by the random walk process. 
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1.01092 
±0.00027 
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±0.00030 


1.01142 
±0.00008 


1.01136 
±0.00010 



Table 2: The thermodynamic hmit of transition temperatures for the — 50 model obtained 
in the methods 'b' and 'c'. Results obtained from scaling of the specific heat peak position 
and fine tuning of free energy till the well depths are equal, are shown in the table along 
with the errors in the linear fits. 
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Figure 1: The histogram fluctuations Aif-^ for the k*^ iteration are plotted against the MC 
sweep index j for the 16 x 16 lattice and — 50. The values of In / are indicated in the 

figures. The histograms are averaged over 50 independent simulations. The error bars are 
shown in the diagram. 

Figure 2: Logarithm of the saturation values of the histogram fluctuation, ln(/3fc) is plotted 
against ln(ln/) for 8x8 and 16 x 16 systems for p'^ = 1. The slopes of the two linear fits 
are given in the text. The error bars are of the dimensions smaller than the symbols used 
for plotting. 

Figure 3: The plots similar to those in fig. 2 forp^ = 50. The error bars are of the dimensions 
smaller than the symbols used for plotting. 

Figure 4: Plot of In A, the logarithm of the MC sweeps necessary to reach saturation of the 
histogram fluctuation 

Figure 5: The change in the value of the logarithm of the density of states, Agk{Ei) that 
occurs during an iteration is plotted against the bin energy The changes are shown for 
four iterations, with the values of In / inscribed in the diagrams. The results are for the 
8x8 lattice and = 50. 

Figure 6: Specific heat per particle obtained by method 'b' plotted against the dimensionless 
temperature for different lattice sizes. 

Figure 7: The free energy A obtained by method 'b' plotted against bin energy per particle 
for all lattice sizes. For clarity only a restricted energy range is shown in the plot. 

Figure 8: Finite size scaling of the peak height of C„ per particle obtained by the methods 
'b' and 'c' and plotted against L^. The error bars are of the size smaller than the symbols 
used for plotting. 

Figure 9: Finite size scaling of the free energy barrier at the transition temperatures by the 
two methods plotted against L. The linear fit is for method 'b' only. The error bars are of 
the size smaller than the symbols used for plotting 



Figure 10: Finite size scaling of the transition temperatures obtained from i) peak position 
of Cy and ii) fine tuning of free energy with temperature plotted against L~^. The straight 
lines joining data points are guide to eye. We have used results obtained from both methods. 
The upper pair of the straight lines are obtained from method 'b' and the lower pair from 
method 'c'. The error bars are of the dimensions smaller than the symbols used for plotting. 



